# Created by Yue Liu
# Modified by Zhixuan Li
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.ticker as ticker

fig=plt.figure();
ax1 = Axes3D(fig);

import numpy as np;
import sys;

list=[];
filename=sys.argv[1];
try:
    fp = open(filename,"r");
    print("%s Open successfully" % filename);
    for line in fp.readlines():
        list.append(line);
    fp.close();
except IOError:
    print("%s Open failed" % filename);

import re;
data = re.split(" ",list[0]);
#for i in range(1,13):
#    data[i] = int(data[i]);

D = int(data[1]);
n = int(data[2]);
q0 = np.array([int(data[3]),int(data[4]),int(data[5])]);
x0 = float(data[6]);
y0 = float(data[7]);
z0 = float(data[8]);
a = float(data[9]);
b = float(data[10]);
c = float(data[11]);
theta = float(data[12]);

ax1.set_xlabel('X');
ax1.set_xlim(0, n);
ax1.set_ylabel('Y');
ax1.set_ylim(0, n);
ax1.set_zlabel('Z');
ax1.set_zlim(0, n);
tick_spacing = 1;
ax1.xaxis.set_major_locator(ticker.MultipleLocator(tick_spacing));
ax1.yaxis.set_major_locator(ticker.MultipleLocator(tick_spacing));
ax1.zaxis.set_major_locator(ticker.MultipleLocator(tick_spacing));

#vmax = np.arccos(min(max(-z0/c,-1),1));
#vmin = np.arccos(min(max((n-z0)/c,-1),1));

PI = np.pi;
u = np.linspace(0, 2*PI, 30);
v = np.linspace(0, PI, 30);

xx = x0 + a * np.outer(np.cos(u), np.sin(v));
yy = y0 + b * np.outer(np.sin(u), np.sin(v));
zz = z0 + c * np.outer(np.ones(np.size(u)), np.cos(v));

if data[0]=='e':
    ax1.plot_surface(xx,yy,zz, cmap='bone', alpha=0.6);
#ax1.contourf(xx,yy,zz,zdir='z', offset=0,cmap="rainbow")  
#ax1.contourf(xx,yy,zz,zdir='x', offset=0,cmap="rainbow") 
#ax1.contourf(xx,yy,zz,zdir='y', offset=0,cmap="rainbow")   

from scipy.special import comb, perm;
ps = [];
pn = comb(n+D,D);
ls = re.split(" ",list[2]);
for i in range(0,len(ls)-1):
    ls[i] = int(ls[i]);

for i in range(0,int(pn)):
    pp = np.array([ls[3*i],ls[3*i+1],ls[3*i+2]]);
    ps.append(pp);

for pp in ps:
    ax1.scatter3D(pp[0],pp[1],pp[2],c='red');

print(len(ps));
ls = len(ps);
for i in range(0,ls):
    for j in range(i+1,ls):
        d = np.linalg.norm(ps[i]-ps[j]);
        if d==1:
            p = ps[i];
            q = ps[j];
            xx = np.linspace(p[0],q[0]);
            yy = np.linspace(p[1],q[1]);
            zz = np.linspace(p[2],q[2]);
            ax1.plot3D(xx,yy,zz,'blue');
            print(i+1,j+1);

e0 = np.array([np.cos(theta)*np.cos(theta) - np.sin(theta)*np.cos(theta)*np.sin(theta), np.sin(theta)*np.cos(theta) + np.cos(theta)*np.cos(theta)*np.sin(theta), np.sin(theta)*np.sin(theta)]);
e2 = np.array([np.sin(theta)*np.sin(theta), -np.cos(theta)*np.sin(theta), np.cos(theta)]);
e1 = np.cross(e2,e0);

print(e0,e1,e2);

ax1.scatter3D(q0[0],q0[1],q0[2], c='black');
xx = np.linspace(x0,e0[0]+x0,20);
yy = np.linspace(y0,e0[1]+y0,20);
zz = np.linspace(z0,e0[2]+z0,20);

if data[0]=='b':
    ax1.plot3D(xx,yy,zz,'red');

xx = np.linspace(x0,e1[0]+x0,20);
yy = np.linspace(y0,e1[1]+y0,20);
zz = np.linspace(z0,e1[2]+z0,20);

if data[0]=='b':
    ax1.plot3D(xx,yy,zz,'green');

xx = np.linspace(x0,e2[0]+x0,20);
yy = np.linspace(y0,e2[1]+y0,20);
zz = np.linspace(z0,e2[2]+z0,20);

if data[0]=='b':
    ax1.plot3D(xx,yy,zz,'yellow');

print(np.dot(e0,e1));
print(np.dot(e1,e2));
print(np.dot(e0,e2));

print(np.arctan2(e1[2],e2[2])*180/PI);
print(np.arctan2(-e0[2],np.sqrt(e1[2]**2+e2[2]**2))*180/PI);
print(np.arctan2(e0[1],e0[0])*180/PI);

plt.show();





             
             

